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cd ! Abstract 

a 

O I In a recent letter [Phys. Rev. Lett. 75, 2360 (1996)] we briefly discussed 

the existence and nature of ferroelectric order in positionally disordered dipo- 
lar materials. Here we report further results and give a complete description 
^vq . of our work. Simulations of randomly frozen and dynamically disordered 



dipolar soft spheres are used to study ferroelectric ordering in non-crystalline 



CN ' systems. We also give a physical interpretation of the simulation results in 

O ' 

[^^ ' terms of short- and long-range interactions. Cases where the dipole moment 



has 1, 2, and 3 components (Ising, XY and XYZ models, respectively) are con- 



C^ \ sidered. It is found that the Ising model displays ferroelectric phases in frozen 

^ \ amorphous systems, while the XY and XYZ models form dipolar glass phases 

^^ , at low temperatures. In the dynamically disordered model the equations of 

O ■ motion are decoupled such that particle translation is completely independent 

O ■ 

of the dipolar forces. These systems spontaneously develop long-range ferro- 

i^i electric order at nonzero temperature despite the absence of any fined-tuned 

^ ' short-range spatial correlations favoring dipolar order. Furthermore, since 

this is a nonequilibrium model we find that the paraelectric to ferroelectric 
transition depends on the particle mass. For the XY and XYZ models, the 
critical temperatures extrapolate to zero as the mass of the particle becomes 
infinite, whereas, for the Ising model the critical temperature is almost inde- 
pendent of mass and coincides with the ferroelectric transition found for the 
randomly frozen system at the same density. Thus in the infinite mass limit 



the results of the frozen amorphous systems are recovered. 
PACS numbers: 64.70.Md, 77.80.-e, 82.20Wt 
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I. INTRODUCTION 

Computer simulations of fluids of strongly interacting dipolar spheres have established 
the existence of a stable ferroelectric liquid crystal phase |jl],|^ . The phase is truly fluid in 
that it exhibits long-range orientational (ferroelectric) order, but only short-range spatial 
correlations. Interestingly, the local spatial correlations found in these ferroelectric fluids 
were similar to those in the ferroelectric tetragonal-I lattice, which is the lattice structure 
believed to be the low-temperature ferroelectric solid phase for dense dipolar hard spheres 
0]. The translational mobility of the particles in the fluid phase allowed specific short-ranged 
correlations to build up, and the system spontaneously polarized. /^From these results, it 
was believed that the establishment of a ferroelectric liquid phase was largely driven by 
well-tuned specific short-range spatial correlations [Q. 

On the other hand, experimental systems of near spherical Fe304 magnetic particles 
in a frozen non- magnetic solvent §-0, which interact via magnetic dipole moments, do 
not exhibit a ferroelectric phase upon cooling, but rather, develop magnetic irreversibilities 
similar to the situation in random magnetic systems known as spin glasses |^. In these 
systems, the particles are frozen at random locations at all temperatures, but above a certain 
temperature the particle dipoles can freely rotate |^ The simplest interpretation of these 
results is that the lack of fined-tuned spatial correlations inhibit the formation of long-range 
ferroelectric order, while the large degree of randomness in the positions lead to random 
frustration for the dipolar interactions and, consequently, to the formation of a "dipolar 
glass" at low temperature. The anisotropic nature of the dipolar interaction would lead one 
to expect that the glass transition in these frozen ferrofluids is of thermodynamic (as opposed 
to purely dynamic) origin, and hence occurs at a nonzero glass transition temperature, Tg, 
characterized by a divergent spin-glass susceptibility @, IC]. 



To summarize, the naive picture that emerges is that some crystalline lattices (e.g., body- 
centered cubic, face-centered cubic, tetragonal-I) and some dipolar fluids can spontaneously 
develop long-range ferroelectric order [|11[] because they have "suitable" spatial correlations. 
On the other hand, low-density randomly frozen dipolar systems lack these correlations. 



exhibit random frustration, and therefore have dipolar glass ground states |T2[. Also, we 
would expect that dilution of the dipoles on a crystalline lattice, which displays ferroelectric 
order for full occupation, will eventually lead to a transition from ferroeletric order to dipolar 
glass below a critical occupation density, again due to the increase build up of random 
frustration with decreasing dipole density P,p!3|. 



However, in recent papers, Zhang and Widom ||14| , ^5| proposed a mean field theory that 
predicts ferroelectric phases in dipolar systems that lacked any specific spatial correlations, 
provided the density of the particles, p, was above a critical value, pc- They considered 
amorphous solids of dipolar hard spheres where the particles were free to rotate, but were 
frozen at random sites. Specifically, they assumed complete randomness where the radial 
distribution function, g{r), describing the probability that two particles are separated by 
a distance r, was set to g{r) = 1 for r > cr, where a is the diameter of the sphere. Their 
prediction of ferroelectric phases in dipolar systems that lack any specific spatial correlations 
suggests that well-tuned short-range structure may not be necessary for ferroelectric phase 
formation. Although the experimental results on frozen ferrofluids seem to contradict this 
assertion, Zhang and Widom suggested that the experimental dipole density was possibly 
too low (i.e., below p^) for ferroelectric order, and hence a dipolar glass state was found. 

The role of short-ranged spatial correlations on ferroelectric phase formation is still not 
well understood. The question of whether or not dense spatially disordered dipolar materials 
can have ferroelectric phases was briefly discussed in Ref. [jl6[. In this paper, the effect of 



spatial disorder on the phase behavior of dense dipolar systems is investigated with molecular 
dynamics (MD) and Monte Carlo (MC) simulations. Systems where the dipole vector has 
one, two and three components are considered and we refer to these as the Ising, XY and 
XYZ models, respectively. The fundamental forces that promote and destroy ferroelectric 
order in dipolar systems are identified and discussed. 

The remainder of this paper is organized as follows. In Section II we briefly discuss the 
generic temperature vs density phase diagram one might expect to find for diluted ferro- 
electric lattices and amorphous dipolar systems. In Section III, the models and simulation 
methods employed are described. Section IV is concerned with ferroelectric and dipolar 



glass ordering found in amorphous frozen dipolar systems. A new simulation technique de- 
vised to study orientational order in spatially random media is introduced in Section V. 
The results of these "dynamically disordered" simulations offer significant insight into the 
observed behavior of both dipolar fluids |l|J^ and the dipolar amorphous solids ||^-|7,1^ 



A mean field theory which describes the ferroelectric transition when only reaction field 
interactions are present is given in Section VI. This helps us understand the competition 
between the long-range reaction field interactions and shorter-ranged contributions to the 
energy. This competition is a key feature determining whether or not a system displays 
ferroelectric order. Finally, our main results and conclusions are summarized in Section VII. 

II. TEMPERATURE-DENSITY PHASE DIAGRAM IN RANDOM DIPOLAR 

SYSTEMS 

It is useful to briefly discuss the temperature, T, vs density, p, phase diagram (sketched 
in Fig. 1) that one might expect to observe in random dipolar systems based on our current 
understanding of randomly frustrated magnetic spin glasses P,p|, p!3| , p!7 . 



First, consider a Bravais lattice (e.g., body-centered cubic) where all sites are occupied 
by an n-component classical dipole such that p = Pmax- We will assume here that the 
positions of the particles at the lattice sites are fixed. At high temperature, the system is 
in a disordered paraelectric phase (see Fig. 1). For a perfect lattice (all sites occupied), a 
transition to a long-range ferroelectric phase occurs at some critical temperature, Tc [ll8| . 
As the system is diluted, a ferroelectric phase remains observable for sufficiently large p > 
Pc but with Tc decreasing as more sites are vacated and p decreases. As well, the zero- 
temperature polarization, P{T = 0), decreases from its full value at p = Pmax as p approaches 



Pc from above. Based on work on spin glasses [£|,|17|; "we expect that upon cooling within the 
ferroelectric phase, there will be a low-temperature "mixed" state {pc < p < Pmax) where 
strong irreversibilities and glassy behavior develops (i.e., the polarization in field-cooled and 
zero field-cooled experiments differs below the long-dashed line in Fig. 1, which is commonly 
referred to as the Almeida- Thouless line in the mean- field theory of Ising spin glasses M). 



However, long-range ferroelectric dipolar order is not lost in the mixed phase and there is 
likely no decrease of the polarization upon cooling from the ferroelectric state down into the 



mixed state ||T^. For p < pc, the random frustration leads to a transition to an unpolarized 
dipolar glass state characterized by an Edwards- Anderson order parameter ^. Most naively, 
we expect that the intrinsic anisotropy of dipolar interactions stabilizes a dipolar glass phase 
at nonzero temperature in the three-dimensional Ising spin glass universality class 0] (with 
the proviso that recent numerical work suggests that the existence of a thermodynamic 
spin-glass transition in the three-dimensional Ising spin glass is still not completely settled 
19|.). In diluted magnetic systems with short-range frustrated interactions, such as diluted 



CdMnTe and EuSrS, one finds that the glass transition temperature vanishes for a nonzero 
density, p_, which corresponds to a site percolation threshold |^. For classical dipoles, p_ is 
zero due to the long-range nature of the interactions. However, for small non-classical dipole 
moments, quantum effects will likely move p_ to a finite value due to the random-transverse 
fields acting on each spin in the frozen state ||20|| . 

Consider now the case of amorphous dipolar systems where there are little or no spatial 
correlations among the locations of the dipoles (i.e., the dipoles are not on lattice sites) 0. 
Recently, Zhang and Widom |1^,|15| have argued using a mean field model that frozen amor- 
phous dipolar systems could display spontaneous ferroelectric order for a packing density 



p > Pc{n) for n = 1 (Ising model) and 3 (XYZ model) component dipoles |jT4|,|T5[. It should 
be noted here that one clearly cannot make arbitrarily dense amorphous systems lacking 
all spatial correlations. However, the critical density for ferroelectric order in amorphous 



Ising and XYZ systems was found in Ref. [14,f5l to be reasonably less than the close-packed 



value, and it should be possible to test the theoretical predictions for a range of densities. In 
other words, varying the temperature for sufficiently large p, one should according to Zhang 
and Widom follow the trajectory a) in Fig. 1, and find a ferroelectric phase below Tc{p,n). 
The main motivation of the present study was to investigate this possibility. Below, 
we report results from extensive computer simulations which show that for a high density 
amorphous system, the Ising model does display a ferroelectric transition. In fact, for the 
density considered, the Ising system developed a polarization close to the full maximum 
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value (i.e., the Ising system considered was at p ^ pc). However for the same and even 
larger densities, we found that the XY and XYZ models showed instead a dipolar glass 
transition and no ferroelectric order was observed. Thus, these systems followed path b) in 
the phase diagram sketched in Fig. 1. 

III. THE MODEL AND SIMULATION DETAILS 

We consider systems of point dipoles embedded at the center of soft spheres. Soft spheres 
(ss) are defined by the pair potential 

nss(12) = 4£((7/r)i\ (1) 

where the parameters e and a are the fundamental units of energy and length, and r is the 
distance between the particles. The dipole-dipole (DD) interaction is given by 

UDDil2) = -3(/Xi ■ r)(/X2 ■ r)/r^ + /Xi ■ Ma/r' , (2) 

where /Xj is the dipole of particle i and r is the interparticle vector. As noted above, we 
consider Ising, XY and XYZ models. In all three cases the potential is given by Eq. (2). 

All calculations were carried out employing periodic boundary conditions and the long- 
range dipolar forces were taken into account using Ewald summation methods. The dielectric 
constant of the surrounding continuum, e', necessary in the Ewald method il| j2l|j2^ was taken 



to be infinity (i.e., conducting boundary conditions). For the present dense strongly dipolar 
systems the dielectric constants are sufficiently large that conducting boundary conditions 
are an appropriate choice. 

The existence of a ferroelectric phase can be detected by calculating the average polar- 
ization per particle, P, defined as 

^=^(EArd), (3) 

where d is a unit vector in the direction of the total instantaneous moment, M = J2f=i A^ij 
and A^ is the number of particles in the system. In ferroelectric systems, P is nonzero and 
tends to one as the ferroelectric order increases. In a disordered phase, (either paraelectric 



or dipolar glass) P will be zero if the system is sufficiently large. However, in simulations P 
must be expected to exhibit significant system size dependence and, as demonstrated below, 
this must be carefully checked. 

In the present work, we are also interested in the possibility that a system could orien- 
tationally "freeze" into an unpolarized state at low temperatures. Such an orientationally 
frozen, or dipolar glass phase can be detected by calculating the root-mean-square (rms) 
dipole length, (S'rms), given by [^ 



(^™s) = ^[E(<«i«>-<™^>)]'^% (4) 

i 

where < rrij >= r^^ J2t'=o AjI'^')) ^^^ r is the number of MD timesteps, or MC sweeps (i.e., 
N attempted moves). Briefiy, for ferroelectric systems both P and (Sj-^s) will be nonzero. 
If P ^ but (S'rms) is nonzero, the system is an orientationally frozen dipolar glass. If both 
P and (S'rms) are near zero, the particles are freely rotating as in a plastic crystal or normal 
isotropic fiuid. 

Systems of dipolar soft spheres can be characterized by specifying the reduced density, 
p* = Na^/V, the reduced temperature, T* = kT/e, where k is the Boltzmann constant, 
and the reduced dipole moment, yU* = {fi^/ea^Y^'^. In the present work, the reduced dipole 
moment and reduced density were constant at /i* = 4 and p* = 0.8. This density is well 
within the range where Zhang and Widom predict a ferroelectric phase. For example, the 
lowest density of the ferroelectric phase is predicted to be p* = 0.31 for the Ising system 
and p* = 0.55 for the XYZ model [0,^. In the present case with /x* = 4 and p* = 0.8, the 
theory of Zhang and Widom predicts that the Ising model will be ferroelectric if T* < 35.2 
and the XYZ model if T* < 4.8 E^ 



IV. RANDOMLY FROZEN SYSTEMS 

In our simulations the "randomly" frozen spatial structure is taken to be a typical fiuid 
configuration of soft spheres at T* = 10.5 and p* = 0.8. Such configurations are readily 
generated by MD simulation of the soft-sphere fiuid. Unfortunately, it is impossible to have 
a truly random and uncorrelated (i.e., the radial distribution function, g{r) = 1 for r ^ cr) 
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spatial configuration at p* = 0.8. The radial distribution function for a soft-sphere fluid at 
p* = 0.8 and T* = 10.5 is shown in Fig. 2. Clearly, the spatial correlations in this system are 
weak and very short-ranged. The point dipoles are located at the centers of the soft spheres, 
initially with random orientations. Constant temperature MD simulations [^] were carried 
out for the XY and XYZ models. The reduced simulation timestep, St* = {e/ma'^y^'^St = 
0.00125, where m is the particle mass, was used together with the reduced moment of inertia 
/* = I/ma'^ = 0.025. Of course, these MD simulations involved only rotational motion since 
the particles remain spatially frozen. Monte Carlo simulations were performed for the Ising 
model. One Monte Carlo "sweep" consisted of A^ attempted random flips of the dipolar 
orientations. Typically, the results reported involved aging runs of 100,000 MD timesteps 
or MC sweeps followed by production runs of the same length. 

The average polarization as a function of 1/A^ for the Ising, XY and XYZ models is 
shown in Fig. 3. The systems considered ranged from 108 to 864 particles. The polarization 
values that are plotted were obtained at the lowest temperature where equilibrium could 
be achieved independent of the starting conflguration. We call this temperature T^i^. The 
values of T^i^^ for the Ising, XY and XYZ models are 10.0, 4.0 and 3.5, respectively. Below 
these temperatures, simulations that were started from perfectly aligned and random states 
did not converge to the same result, at least not in simulation runs of practical length. 
However, the values of T^^j^ attained are within the temperature range where Zhang and 
Widom predict a ferroelectric phase. As discussed below, it appears that the Ising system 
develops ferroelectric order with P > in the thermodynamic limit. In this case, T^^^ may 
qualitatively correspond to the long-dashed line in Fig. 1. 

/^From Fig. 3 we see that the Ising model at T^^j^ = 10.0 is almost completely polarized 
and shows little or no system size dependence. A detailed P vs T* plot for the frozen A^ = 256 
Ising system is given below (see Fig. 7). It can be seen that ferroelectric order develops 
spontaneously in this system at T* ^ 25. This transition temperature is signiflcantly lower 
then that (i.e., T* = 35.2) predicted by the theory of Zhang and Widom. Returning to Fig. 
3, we see that the XY and XYZ models at T^^^ show signiflcant polarization for the 108 
and 256 particle systems. However, in both cases the polarization decreases monotonously 



with increasing system size and appears to approach zero in the thermodynamic hmit. The 
observed polarization of the XY model is strongly dependent on system size, decreasing 
from ~ 0.49 for A^ = 256 to ~ 0.10 for A^ = 864. The polarization for XYZ model shows a 
similar, although not as pronounced, system size dependence. 

Further information about the behavior of the XYZ and XY models can be obtained by 
examining the temperature dependence of (5'rms)- Results for the XYZ model (A^ = 256) 
are shown in Fig. 4. These results were obtained for samples initially begun with random 
dipolar orientations. We see that (^rms) is zero at high temperatures as expected, but 
becomes nonzero and grows with decreasing T* at lower temperatures. Similar behavior 
was observed for the XYZ system with 864 particles, and for the XY model. The growth of 
('S'rms) at low temperatures without a corresponding development of polarization provides 
qualitative evidence that the XYZ and XY models freeze orientationally to form dipolar 
glasses. 

To summarize, we find no evidence of ferroelectric states in the thermodynamic limit for 
the randomly frozen XYZ and XY models. This clearly disagrees with the theory of Zhang 
and Widom, which predicts that the XYZ model should have a stable ferroelectric phase in 
the temperature range we consider. Additional MD calculations at p* = 1.05 were carried 
out for the XYZ model, but again no ferroelectric behavior was observed. The Ising model 
does have a ferroelectric phase, however, spontaneous polarization was observed at T* ^ 25, 
which is much lower than the transition temperature predicted by Zhang and Widom. 

V. DYNAMICALLY DISORDERED SYSTEMS 

In the simulations described above, a configuration was selected from a MD simulation 
of soft spheres at p* = 0.8 and T* = 10.5, and was then used as a "typical" randomly frozen 
system. Point dipoles were embedded at the centers of the soft spheres and rotational MD 
or MC simulations were performed. Ideally, as in spin glass models P], many (i.e., 100-1000) 
randomly frozen configurations should be used to obtain accurate values of the "disorder- 
averaged" polarization. However, for the dipolar systems considered here this would be a 
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very laborious procedure requiring many long simulations. 

As an alternative approach to the study of orientational ordering in random media, we 
have used a MD simulation technique where the rotational and translational equations of 
motion are completely decoupled. That is, we consider dipoles embedded in a dynamic 
random "substrate" rather than in a frozen system. The underlying soft-sphere substrate is 
a simple fluid. It has no long-range positional correlations and the short-range correlations 
are not influenced by the dipolar interactions (i.e., the spatial correlations are identical to 
those of a soft-sphere model [Eq. (1)] at a given density and translational temperature). 
However, the "equilibrium" state of the dipoles will depend on the underlying motion of 
the substrate. This model is similar in spirit (but not equivalent to due to its lack of 
obvious energy currents) to those used in recent studies of non-equilibrium phase transitions 



in magnetic systems subject to Levy flights ||26|. It is fundamentally different from dipolar 



fluid simulations in that the spatial structure is not affected by dipole-dipole interactions, 
but it is also not an amorphous solid simulation since the particles move. 

With this technique, we can gain a good deal of insight into the role of short-ranged 
spatial correlations on phase behavior. By controlling the rate that the substrate moves 
relative to the rate of dipolar reorientations, we can effectively "turn on" or "turn off" the 
specific details of the random spatial structure which the dipoles "see" . The implementation 
is straightforward. The force between two particles is simply given by 

f(12) = -VrMss(12) , (5) 

where Mss(12) is the soft-sphere potential defined above. The torques are given by the 
dipole-dipole interactions. The spatial structure of the system is then determined only by 
the soft-sphere part of the pair potential. The translational and rotational temperatures are 
decoupled such that the structure and motion of the substrate is not affected by changes in 
the rotational temperature. The rate at which the substrate moves relative to the dipolar 
reorientations can be varied by changing the particle mass. Obviously, since it is a classical 
fluid, adjusting the particle mass will have no effect on the equilibrium structure of the 
soft-sphere substrate. For large masses, the substrate changes slowly relative to dipole reori- 
entations. An extrapolation to the inflnite mass case would give results of a randomly frozen 

11 



system. For light masses, the substrate moves rapidly relative to dipole reorientations. The 
substrate motion may be so rapid that the dipoles cannot react to structural details and a 
mean-field- like limit is reached p6|. Thus, from a dipole's point of view, increasing the mass 
effectively "turns on" the specific structural details of its surrounding dipolar environment. 
The moving substrate is a means of simulating dipolar systems in a dynamically random 
medium that lacks any specific spatial correlations. Clearly, in these systems positional 
correlations which favor ferroelectric order are not present initially and, due to the decou- 
pling, such correlations cannot develop as the system evolves. Contact with the randomly 
frozen systems can be made by examining the behavior at intermediate masses and then 
extrapolating to the infinite mass limit. 

In the decoupled MD simulations, the underlying substrate is a soft-sphere fiuid at p* = 
0.8 and T*(translational) = 10.5. All decoupled simulations were carried out using the 
reduced timestep, St* = {e/m'a'^Y^'^5t = 0.00125, and the reduced moment of inertia /* = 
I/m'a'^ = 0.025. The equations of motion of the soft-sphere substrate were written in terms 
of the reduced mass, m* = m/m', which could be varied to change the rate of translational 
motion of the substrate. Note that I* and the rotational equations of motion which govern 
the dipoles have no dependence on m*. 

P versus T* (rotational) for the XYZ model is plotted in Fig. 5. Systems with particle 
masses m* = 1, 5 and 10 are included. Clearly, spontaneous polarization develops for all 
systems and the temperature at which P begins to rise decreases with increasing mass. For 
m* = 5, the transition occurs at T* ^ 0.55, and for m* = 10, at T* ^ 0.25. For m* = 1, we 
have also calculated the Binder ratio, SfBinder, defined as ^ 

^Binder = (5/2) - {3/2){\M\')/{\M\Y, (6) 

for systems with 64, 108 and 256 particles. A plot of f^Binder vs T* is included as an inset in 
Fig. 5. A clear crossing, and hence a thermodynamic transition, at T* ^ 1.9 is evident. 

XYZ systems with larger masses (up to m* = 20) were investigated but all were dis- 
ordered above T* = 0.1. Below T* = 0.1 calculations for very large masses converged too 
slowly to be useful. Similar results were obtained with A^ = 108, 256, and 500 particle 
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systems and the polarization showed no significant dependence on system size. The results 
shown in Fig. 5 strongly suggest that for any finite mass the XYZ model will spontaneously 
polarize at some rotational temperature, but as the mass becomes very large the transition 
temperature will approach zero. As shown in Fig. 6, the dynamically decoupled XY model 
behaves much as the XYZ system. It can be seen that systems with m* = 1, 5 and 10 
spontaneously polarize at T* (rotational) ~ 6, 2, and 1.8. Again, the transition temperature 
decreases with increasing mass. 

In the Ising model the potential does not vary as a continuous function of orientation 
and hence this model is not well suited to MD simulations. Therefore, a MC scheme was 
devised that allowed the substrate to move independently of the Ising dipoles. This in- 
volved combining a soft-sphere MD simulation with a MC Ising dipole simulation. The 
soft-sphere substrate evolved as in the XY and XYZ decoupled calculations, however, after 
each MD timestep, A^ attempted dipole fiips (one MC sweep) were performed using the usual 
Metropolis Monte Carlo method p5[. In Fig. 7, P vs T* (rotational) results are plotted for 
m* = 1, m* = 5 and the randomly frozen system (i.e., m* = cxd). We see that the ordering 



behavior of the Ising model is essentially independent of mass [^] and that the results for 
the dynamically disordered systems lie very close to those for the randomly frozen case. 
Reduced constant volume heat capacities per particle, Cy/Nk, obtained by numerically dif- 
ferentiating the average dipolar energy with respect to the rotational temperature are also 
shown in Fig. 7 (see inset) |2^. The randomly frozen and m* = 5 results are very similar 



and indicate a phase transition at T* ^ 25. 

The dependence of the ferroelectric transition temperature, Tp (rotational), on particle 
mass m* for the XY and XYZ models is shown in Fig. 8. The transition temperatures were 
estimated from heat capacities (see Fig. 8, inset) obtained from numerical differentiation 
of the dipolar energy per particle. Results for the Ising system are not plotted because the 
transition temperature is essentially independent of the mass [0. As the mass increases 



the transition temperature drops for both the XY and XYZ models. As noted earlier, for 
large masses and low rotational temperatures convergence becomes prohibitively slow, but 
it seems reasonable to assume that the graph would simply continue with the transition 
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temperature approaching zero in the infinite mass hmit. This is clearly consistent with 
the fact that we did not find a ferroelectric phase for randomly frozen systems at finite 
temperatures. 

VI. MEAN FIELD THEORY: ISOLATING THE EFFECT OF LONG-RANGE 

INTERACTIONS 

In order to better understand the simulation results discussed above, it is useful to con- 
sider a simple mean field theory for a system where only long-range reaction field interactions 
are included. The properties of such a system are completely independent of spatial struc- 
ture and it is instructive to compare its behavior with the simulation results for the various 
models. 

We consider N particles in a spherical cavity of radius, Tc, surrounded by a continuum 
of dielectric constant e . The reaction field, R, within the cavity arises from the polarization 
of the surroundings by the dipoles in the sphere and is given by pSJpQ 



R = fie')M/rl , (7a) 

where 

M = J2fi,, (7b) 

j=i 

is the total dipole moment of the sphere and /(e') = 2(e' — l)/(2e' + 1). If, as in the 

simulations, e' = oo, then /(e') = 1 and we shall write all further expressions for this 

specific case. It is also obvious from the definition of /(e'), that as e' increases /(e') rapidly 

approaches 1 and becomes effectively independent of the exact value of e'. This is the 

justification mentioned above for using e' = cxd in simulations of dipolar systems which large 

dielectric constants. 

The total instantaneous energy of the system, U, can be written in the form 

1 ^ 
[/ = --Em.-R, (8) 
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where factor of 1/2 arises when one calculates the energy of a point dipole in a reaction field 
P3| , |2D|J5(]| . It takes account of the work required to polarize the surrounding continuum. 
Clearly, the low energy state is when the dipoles are aligned with the reaction field. Using 
Eqs (7), U can be expressed in the form 

u = -^p^^ J^2^2^^^i■^^j , (9) 

where we have also used N = iirr^p/S and p = N/V. 

In the canonical ensemble the average polarization can be obtained by evaluating 

where Pins = (1/^) J2iLi Aj ' d is the instantaneous polarization and dfl^ represents in- 
tegration over the angular coordinates of the N particles (here d defines the z axis of the 



coordinate system) . If we make the mean field approximation |^ | 



(A. - Pd) ■ (A, - Pd) = , (11) 



then Eq. (9) simplifies to 



U^-^pp'pJ:^rd+^p'NP\ (12) 

and Eq. (10) can be easily solved for P. Actually, here we dealing with an effective infinite 
range system where each dipole is coupled with every other dipole with the coupling constant 
27Tpp'^/3N [see Eq. (9)]. In this infinite range case, the mean field solution is equivalent to 
the Hubbard-Stratonovich solution and is essentially exact [^. For the XYZ, XY and Ising 
models, respectively, one obtains 

P = coth(x) - 1/x , (13a) 



Io(a;) 



P = tanh(a;) , (13c) 
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where x = 4npf3^'^P/3 and Iri(x) is the modified Bessel function of order n. The average 
energy per particle is given by 

{U)/N = -2npi2^Py3 , (14) 

and the constant volume heat capacity can be obtained by taking the differentiation with 
respect to temperature. 

The mean field theory predicts ferroelectric transitions for all three models. By expanding 
Eqs (13) about x = one finds that the mean field ferroelectric transition temperatures are 
given by 

T; = 47rp*fi*^/3n , (15) 

where n is the number of dipole components. For the Ising and XYZ models Eq. (15) agrees 
with the "conventional mean field" results given by Zhang and Widom [|T3|. For the present 
parameters {p* = 4, p* = 0.8) the transition temperatures obtained are 17.9, 26.8 and 
53.6 for the XYZ, XY and Ising models, respectively. We emphasize that the driving force 
for the transition in this simple theory is just the reaction field due to the self-consistent 
polarization of the surrounding continuum. 

It is clear from the above discussion that the reaction field interactions favor ferroelectric 
order for all three models. However, as described in Sections IV and V, computer simulations 
of randomly frozen or dynamically decoupled systems indicate the existence of ferroelectric 
order only for the Ising model. Ferroelectric phases were not observed for the XY and XYZ 
models and the MD evidence strongly suggests that these systems are unpolarized at all 
finite temperatures. This behavior can be explained if we consider that the reaction field, 
R, is only one contribution to the total local field, Eiocai, experienced by a particle. The 
local field can be written as 

Elocal = R + E , (16) 

where the remaining contribution, E, is dependent on positional correlations and may or 
may not favor ferroelectric order. If R dominates, ferroelectric phases are to be expected. 
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However, if the local field is largely determined by E, then the existence, or non-existence of 
ferroelectric phases will depend on the details of the spatial correlations (eg., as in Bravais 
lattices). 

/,From this point of view, it is useful to divide the average energy obtained in the simu- 
lations into reaction field (RF) and structurally dependent (SD) parts such that 

{U) = (f/Rp) + (f/sD) . (17) 

Of course, (f/Rp) depends only on particle orientation with respect to R, whereas {Uq,d) will 
be sensitive to details of the local spatial structure experienced by a particle. The total 
energy and both contributions obtained in dynamically decoupled simulations of all three 
models are shown as functions of T* (rotational) in Fig. 9. Results for different values of 
the reduced mass are included. For all three models we see that (f/Rp) dominates in the 
ferroelectric state and that (f/so) dominates in the paraelectric phase. In fact, in all cases, 
as the system becomes ferroelectric (Usd) increases (i.e., becomes less negative) indicating 
that the structurally dependent interactions do not favor ferroelectric order. 

For the XYZ (Fig. 9a) and XY (Fig. 9b) models the different contributions to the energy 
are strongly mass dependent. As the mass increases, and the particles move more and more 
slowly, (Usb) dominates at lower and lower rotational temperatures. Thus, if we begin at 
a fixed rotational temperature in the ferroelectric phase and increase the mass, the dipoles 
eventually respond to the details of the random local structure, (Usd) dominates, and the 
ferroelectric phase disorders. The Ising model does not exhibit this behavior. In the Ising 
case both contributions to the energy are largely independent of mass (Fig. 9c) and a stable, 
mass-independent, ferroelectric phase is observed. The difference in behavior of the Ising 
model must arise from the fact that, with only two dipolar directions available, its response 
to local structure is much more limited than that possible for the XY and XYZ systems. 

VII. SUMMARY AND CONCLUSIONS 

In this paper, we have used computer simulation methods to explore the phase behavior 
of spatially disordered dipolar systems. Both randomly frozen and dynamically disordered 

17 



models were considered. It was found that the behavior observed depended upon the num- 
ber of components included in the dipoles. Ferroelectric phases were not observed for the 
randomly frozen XY and XYZ models. Rather, these systems appear to freeze orientation- 
ally into dipolar glasses as the temperature is lowered. The behavior of the randomly frozen 
Ising model was significantly different. This model does have a stable ferroelectric phase at 
sufficiently low temperatures. 

We also considered dynamically disordered systems simulated such that the translational 
motion of the particles was independent of the dipolar forces. The rate of translational 
motion was adjusted by varying the particle mass. Interacting dipoles embedded in this 
moving substrate can respond to the short-range spatial structure only if the motion is 
sufficiently slow. For small masses (rapid translational motion), ferroelectric phases were 
observed for all three models. However, as the mass is increased and the dipoles become 
"aware" of the short-range random structure, the ferroelectric order is destroyed in the XY 
and XYZ models. In the infinite mass limit both the XY and XYZ models appear to be 
paraelectric at all temperatures. This is consistent with our observations for the randomly 
frozen systems. In contrast, the Ising model exhibits a paraelectric-to-ferroelectric transition 
at a temperature which is essentially independent of mass and agrees with the result for the 
randomly frozen case. 

In order to understand our observations, it is useful to divide the total local field expe- 
rienced by a dipole into two parts. These are, a reaction field contribution which has no 
dependence on the local spatial structure, and a structure-dependent part which includes 
everything else. Clearly, the average total energy can be divided into corresponding reaction 
field and structure-dependent contributions. Both contributions to the total energy were 
evaluated in our simulations. In all systems where ferroelectric order exists, it was shown to 
be stablized by the reaction field interactions. In spatially random systems, the structure- 
dependent contribution never favors ferroelectric order. In the XY and XYZ models, the 
structure-dependent part dominates and paraelectric dipolar glasses rather than ferroelec- 
tric phases are observed. For the Ising model, the reaction field contribution dominates and 
ferroelectric order is observed. In all likelihood, the Ising model differs from the XY and 



XYZ systems simply because with a one-component dipole the opportunities for strong local 
interactions are severely limited. 
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FIGURES 
FIG. 1. A sketch of a possible phase diagram for spatiahy disordered dipolar systems. The 

various terms, symbols and lines are referred to in the text. 

FIG. 2. The radial distribution function, g{r), for soft spheres at p = 0.8 and T* = 10.5. 

FIG. 3. The polarization P at T*^:^^ vs 100/A^ for the randomly frozen Ising (circles), XY 
(squares) and XYZ (triangles) models. Results are included for 108 (XYZ only), 256, 500 and 864 
particles. 

FIG. 4. The root-mean-square dipole length, (S'rms)) for the frozen XYZ model with N = 256. 

FIG. 5. P vs T* (rotational) for dynamically random XYZ systems. The squares, triangles 
and circles are for m* = 1,5 and 10, respectively. The error bars represent one estimated standard 
deviation. (7Bindcr vs T* (rotational) is shown in the inset for A^ = 64 (squares), 108 (triangles) and 
256 (circles) particles. 

FIG. 6. P vs T* (rotational) for dynamically random XY systems. The squares, triangles and 
circles are for m* = 1, 5 and 10, respectively. The error bars represent one estimated standard 
deviation. 

FIG. 7. P vs T* (rotational) for the Ising model. Results are shown for dynamically random 
systems with m* = 1 (squares) and m* = 5 (triangles) and for the randomly frozen case (solid 
circles). The reduced heat capacities per particle, Cy/Nk, are plotted vs T* (rotational) in the 
inset. 
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FIG. 8. The mass dependence of the disordered-to-ferroelectric transition temperature Tp 
(rotational). The squares and triangles are results for the XY and XYZ models, respectively. The 
values of Tp were obtained from plots of the heat capacity, Cy/Nk, vs T* (rotational) and a 
typical example is shown in the inset. The error bars represent estimated uncertainties in the peak 
position. 

FIG. 9. The contributions to the average energy for the XYZ (a), XY (b) and Ising (c) 
models. T* is the rotational temperature. The solid, crossed and open symbols denote (U), (f/so) 
and (t/Rp), respectively. The squares, triangles and circles indicate m* = 1, 5 and 10, respectively. 
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